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Abstract 



A method for performing high order mesh refinement multigrid computations is pre- 
sented. The Full Approximation Scheme (FAS) multigrid technique is utilized for a sequence 
of nested patches of increasing resolution. Conservation forms are generated on coarse scales 
by additional defect correction terms which counter the local excess fluxes at the boundaries. 
Formulas are given for arbitrary order, extending the existing technique of Bai and Brandt. 
Test calculations are presented for a singular source in three dimensions which illustrate the 
multigrid convergence properties, numerical accuracy, and correct order of the approach. 
Applications for all electron quantum chemical computations are discussed. 
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1 Introduction 



Many if not most computational physics and chemistry problems require consideration of 
a large range of length scales. In protein folding, both short and long ranged interactions 
and competitions between them lead to the final configuration of the molecule. |], [| For 
the interaction between a protein and a nucleic acid or a charged interface, specific ionized 
groups may contribute significantly to the binding, while distant regions of the molecule 
have less importance. || If one studies the electronic structure of large molecules, there is a 
concentration of electron density around the atomic nuclei and between atoms in the chemical 
bonding regions, while often large portions of space exhibit very low and smoothly varying 
density. |3J In numerical simulations of fluid dynamics, there may exist specified regions which 
require a higher resolution treatment locally. || 

Multiscale methods provide one approach for tackling computational problems exhibit- 
ing a range of length scales. || ^] These methods were developed in order to overcome conver- 
gence problems in iterative solutions to partial differential equations. By utilizing approx- 
imations from coarser grids, components of the error on a wide range of length scales can 
be decimated, typically leading to linear scaling computing time with system size. The un- 
derlying differential equations can be represented in various ways including finite differences 
and finite elements. For the present method, high order finite difference representations are 
employed. 

Based on the physical examples given above and many others requiring variable resolu- 
tion, it is appropriate to develop mesh refinement strategies in which fine gridding is focused 
only in those spatial regions which require it.|| The goal is to maintain the linear scaling 
property of the multigrid method while minimizing the prefactor in the scaling relation. The 
computational philosophy adopted here is to generate a sequence of nested regular grids; 
this strategy allows one to use the existing multigrid routines over the mesh refinement 
patches with no significant changes. It also allows for the implementation of accurate high 
order difference equations on a composite mesh. This approach is to be contrasted with the 
many curved grid generation techniques [1C] which have been widely applied in engineering 
applications. 

Bai and Brandt || addressed several pivotal issues concerning extension of the FAS 
multigrid technique for locally refined meshes. First, they developed a A-FMG altorithm in 
order to restore linear scaling behavior which can be lost when many levels of refinements are 
used and thus the coarser global grids are themselves visited in a way which scales linearly 
with the number of levels. The work to accuracy exchange rate A is the Lagrange multiplier 
for the grid optimization equations. Second, they determined that the same interpolation 
order can be used at the mesh boundaries as is used over the rest of the domain. Third, 
they found that local relaxation sweeps near structural singularities (to be differentiated 
from source singularities) can restore convergence rates to those observed away from the 
singularities. Finally, a second order conservative differencing method was developed for 
interior source singularities, when it was found that the only important factor for obtaining 
accurate solutions far from the singularity was to correctly reproduce the source strength 
around that singularity. These mesh refinement techniques were tested on several model two 
dimensional problems. 
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In the present work, the Bai and Brandt conservative differencing method is extended 
to high orders and three dimensions. The method is then successfully tested on a three 
dimensional source singularity, the Coulomb potential. First, a review of the FAS multigrid 
method is given. Next, a simple approach for generating the high order formulas for the 
Laplacian and the interpolation near the boundaries is reviewed. Then, the conservative 
FAS differencing forms are derived from general considerations of balancing the fluxes at 
the boundaries. The three dimensional conservative forms are obtained simply by locally 
averaging the one dimensional fluxes over surfaces at the boundaries. With the inclusion 
of the high order flux corrections, the summation of the defect correction over the mesh 
refinement is zero to machine precision. The boundary corrections lead to accuracy within 
the patch which is the same as that for a high order uniform fine mesh covering the whole 
domain; with no corrections, serious errors occur over the whole composite domain. The 
accuracy outside the refinement zone is improved over that for the high order uniform coarse 
mesh over the whole domain. The correct high order behavior is thus obtained over the 
entire composite domain. The only exceptions are the points one grid spacing away from 
the singularity, where high order is not obtained on uniform nor composite meshes. 

High orders are often required in order to obtain numerically accurate results on three 
dimensional grids of reasonable size. While the methods derived here are general for elliptic 
problems, part of our motivation stems from the recent development, along with other groups, 
of ab initio multigrid methods for quantum chemistry, fill , HH, |13| In our all electron approach, 
all particles are represented numerically on the grid, including the nuclei. In previous work, 
we have carried out density functional calculations on atoms and small molecules on uniform 
grids, and it is clear that the majority of the numerical errors originate from the regions 
around the nuclei. These errors are due both to the finite size of the nucleus on the grid and 
to the poor numerical representation of the core electron orbitals. The methods developed in 
this paper will be incorporated in our electronic structure codes both for the Poisson solver 
and for the eigenvalue solver. 



2 Full Approximation Scheme Algorithm 

The FAS multigrid technique allows for solution of nonlinear problems and is ideally suited 
for the mesh refinement methods presented below. || |8|, || A two dimensional schematic of 
the composite mesh examined here is presented in Figure 1, with two nested patches within a 
full domain. In the test calculations here, four or five levels were employed with the coarsest 
three covering the entire domain (see Figure 2). The Poisson equation to be solved on the 
finest patch of the four level problem is written as: 

jji&jjhA = fhi ^ 

For this case, L M is the finite difference Laplacian on the finest scale, U M is the exact grid 
solution for the potential on that scale, and / is — 4n times the charge density for a three 
dimensional problem. The current approximation to the exact solution will be written in 
lower case: for example the approximate M solution is u . On the next coarser level h3, 
the level M patch covers only a portion of the level h3 domain. The equation to be solved 
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on the h3 level is: 



L h3 U h3 = lfj M + r h3 , (2) 

where is the full restriction operator which performs a weighted local average of the fine 
scale function, and r h3 is the level h3 defect correction given by: 

T h3 = L h3jh3 u M _ jh3 L M u M (3) 

The fine scale function can then be corrected as follows: 

and further iterations are subsequently performed on the hA level. Here I^f is the interpo- 
lation operator. One way to understand the function of the defect correction is to observe 
that, if the exact solution U hA were passed to level h3, no correction would be made. That is, 
the defect correction causes the level h3 equation to 'optimally mimic' the level h4 problem. 
Note that r h3 is only defined over the coarse grid points within the interior region of the H4 
level patch, with zero values outside. 

When using multiple scales, the defect correction includes an additional contribution 
from the previous scale. Here the example is given of the level h2 defect correction computed 
during the final level h4 V cycle: 

T h2 = L h2jh2 u h3 _ jh2 L h3 u h3 + j^hS (5) 

By performing these coarse grid correction cycles recursively to coarser and coarser scales, 
errors of all wavelengths can be effectively removed with only several iterations necessary on 
the fine scale. 

Before deriving the high order conservative forms for the composite mesh computations, 
a simple general procedure for obtaining the required high order Laplacian and interpolation 
operators is summarized. 

3 Direct Method for High Order Formulas 

For completeness, a summary is given here of a general method for developing any of the 



formulas required in a high order multigrid method. Hamming||14|| outlines a direct method 
for obtaining numerical formulas of arbitrary order. Sample a function at N points. Define 
the Lagrange sample polynomials: 

TTi(x) = (x — xi)(x — x 2 ) ■ ■ ■ (x — Xi-i)(x - X i+ i) ■ ■ ■ (x — Xjv), (6) 
which are of order N — 1 in x. Then 

7T J (X,)=0, (7) 

and 
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TTjOi) ^ 0. 

A polynomial which passes through the N sample points yi is: 

7Ti(x) 



N 



P N -l(x) =J2Vi 

i 

Expand the sample polynomial as follows: 



(9) 



JV-1 



TCi{X) = Ci 5 i + Ci^X H hC^ArX 1 ' \ (10) 

where c^n = 1. Then it is easy to see that the product of the coefficient matrix and the 
Vandermonde matrix of the sample points (second matrix below) is the matrix [^(xj)]: 
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7T; 



Thus, except for a normalization factor of 7Tj(xj) for each row of the coefficient matrix, the 
left hand matrix is the inverse of the Vandermonde matrix. 

In the direct method, Hamming shows that there is a simple connection between the 
Vandermonde matrix, the desired weights for the approximation, and a vector of 'moments': 



Xw = m 



(12) 



These moments result from allowing the operation of interest to act on the sequence of 
functions: 1 . For example, for the second derivative operator centered at 

x = 0, the first four elements of the moment vector are 0, 0, 2, 0. Similarly, moments can 
also be obtained for the operations of integration and interpolation. For integration, the 
moments are the integrals of l,x, x 2 , . . . over the sampling domain, and for interpolation, 
they are these elementary polynomials themselves. 

Since the normalized version of the coefficient matrix in Eqn. 11 is the inverse of 
X, the weights for each of the approximations in the multigrid process can be calculated 
to any desired order by one matrix-vector multiply. The normalized coefficient matrices 
for N sampling points are obtained by expanding the sampling polynomial (Eqn. 10) and 
dividing each row by the normalizing factor 7Tj(xj). They are termed 'universal matrices' 
by Hamming due to their generality, i.e. they depend only on the sampling points, not 
on the formula to be approximated. The matrices are tabulated in Ref. |14[] up to seven 
sampling points, which allows for computation up to the 6th order Laplacian. In the present 
work, simple C codes have been written which generate the eight, nine, etc. point matrices 
as well. The weight vectors for the Laplacians through 8th order are given in Table 1. 
The three dimensional versions are generated from the sum of the three orthogonal x, y, z 
axes. We have utilized Laplacians up to 8th order in previous multigrid work on uniform 



domains. [|l^] The interpolation weight vectors for even numbers of sampling points are listed 
in Table 2 through 8th order. Similary, formulas for any higher orders can be obtained 
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from the universal matrices. The high order interpolation formulas only need be used when 
setting the function values on and outside the refinement boundaries which are fixed during 
iterations over the patch. It is essential that the order of the interpolation match the order 
of the Laplacian at these boundaries. Lower order interpolations are adequate during the 
rest of the multigrid processing. 



4 High Order Conservative Forms 

When solving for the potential on coarser scales which contain a mesh refinement patch at 
the next finer level, it is clear from Eqn. 2 that, if the sum of r over the interior domain is not 
zero, additional sources have been introduced. This is in fact the case, which can be shown 
by examination of the r terms in a one dimensional example; most of the interior terms do 
cancel, but nonzero contributions remain at the patch boundaries. The terms which remain 
are of the form of one dimensional flux operators. Without correcting for these new sources, 
the solution will be polluted over the whole domain. The method of Bai and Brandt corrects 
for these sources by introducing local opposing fluxes at the boundary. In this section, their 
second order method is extended to high orders. 

First, the problem is illustrated schematically by using a continuous notation (in the 
grid notation, all integrals go over to sums). The coarse scale is labelled by H and the fine 
by h. It is desired to satisfy: 



Id 



T H dV = 0. (13) 

ID 

where the integration is over the whole patch domain D, including the boundaries. However, 
it is true that 



l^dV^O. (14) 
This integration is only over the interior region of the patch. Therefore: 

f r H dV = [ r t H nt dV + / r b H dV = 0. (15) 
Jd Ji Js 

Here r b H is a boundary term designed to oppose locally the additional terms due to non- 
conservation of source and the S integration is over a narrow strip at the surface. This 
implies: 

J s r b H dV = -J i r t H nt dV. (16) 

The form of T? t is the difference of the Laplacian acting on the coarse scale function minus a 
local average of the Laplacian acting on the fine scale function, Eqn. 3. Therefore, converting 
a volume integral into a surface integral: 

[ r h H dV = - f [(V 2 ) V - ((V 2 ) h u h )]dV = - f [V fe V - V h b u h ]da. (17) 
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The brackets () signify a local average (restriction) of the fine scale Laplacian acting on the 
function, and the gradient operators are obtained by noncancellation of terms near the 
boundary of the volume integral. The final expression shows that the boundary r b H generates 
a flux which locally opposes the flux from the additional sources in the interior. Therefore, 
after collecting the correct units from the two scales, it is apparent that the form for r b H is: 

_ r b H H 2 a = [V^u H - V h b u% (18) 

where H is the coarse grid spacing, a is the numerical prefactor to the Laplacian (see Table 
3), and now the gradients are one dimensional operators directed outward from the surface 
(determined below). Here the gradients simply represent the unitless coefficients since H 2 
and a have been moved to the other side. For example, on a one dimensional domain 
(corresponding to the 2nd order Laplacian): 

[V H u% =uf- uf +1 (19) 

on the left boundary, and 

[V H u% = uf - uf_, (20) 

on the right. 

Since the process of full restriction is a weighted local average over the fine scale function 
(27 points in three dimensions), the averaging in Eqn. 17 can be viewed as follows. First, 
average over the direction normal to the boundary surface, compute the two gradient terms 
on the rhs of Eqn. 18, and then average over the other two directions. The full restriction 
weights for this process along one fine scale dimension are w = [1/4 1/2 1/4]. There is 
no requirement for high order restriction operators, as long as the same restriction method 
is used consistently. Therefore, the coefficients for the two gradient terms in Eqn. 18 can 
be determined by solving the one dimensional problem. In two dimensions, the fine scale 
gradient operator is averaged over three points along the boundary line with weights [1/4 
1/2 1/4]. This yields a local average of the fine scale flux through the boundary. A similar 
procedure was followed in the work of Bai and Brandt. || In three dimensions, the local flux 
average is over a square centered on the location of the coarse scale gradient. The weights 
are 1/4 for the center, 1/8 for the edges, and 1/16 for the corners. 

The one dimensional version of the flux difference in Eqn. 18 was solved for high orders 
by examination of the cancellation of terms near the boundary. The result for the left hand 
side of the coarse scale gradient on a left boundary is given by: 

i 

d-n L +i = c -n L +j i = 0,n L -l, (21) 

3=0 

where til is the number of points in the Laplacian to the left of the center and the C- nL+ j 
are the Laplacian coefficients from Table 1. The rhs side of the gradient is antisymmetric 
with respect to these coefficients. For a right side boundary, all the signs are reversed. 
The locally averaged (in one dimension) fine scale gradient coefficients are: 
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i-1 

e -n L +i = ^ 2c_„ i+ j + C_ nL+ j Z = 0, Ul — 1. (22) 

i=0 

For the fine scale coefficients, the central term always cancels completely, so both gradient 
operators and are centered about the fine grid location one point inside the patch 
boundary All of the coefficients up through 8th order are listed in Table 3, and the terms 
for a 6th order left boundary are shown in Fig. 3 to illustrate the locations. Similarly, 
conservative forms can be derived for higher orders if desired. 



5 Computational Details and Numerical Results 

The computational test case presented here is for the 4th order form. The three coarsest 
scales covered the whole domain, while the finest one or two were nested patches. On the 
three coarsest (full domain) scales, the boundaries were set by fixing the potential at the 
analytical value for a singular source in three dimensions <p(r) = l/r. The boundary was 
fixed with one additional term outside the physical boundary since the Laplacian has two 
terms beyond the center in one dimension. Iterations were performed over all the interior 
points of the full domain or patch. The FAS-MG technique was used in the form of the series 
of nested V cycles as shown in Fig. 2. SOR iterations were employed for all relaxation steps, 
with uj = 1.2. The optimal relaxation parameter was determined empirically for the high 
order case. Full weighting restriction and linear interpolation were used, except 4th order 
interpolation was performed over the patch regions, including the required points beyond 
the boundaries. These points were set such that the Laplacian and defect correction were 
defined over the entire interior of the patch. The boundary potential terms for the patches 
were reset during the correction step after each visit to coarser scales. 

The code was written in C with double precision arithmetic, utilizing the prescription 
of Ref. [HJ for dynamic memory allocation. The test calculations were run on a Pentium 133 



MHz processor laptop with 40Mb of RAM, requiring a total of roughly 15 relaxation sweeps 
on the fine scale and 3 seconds of total processor time for convergence. The 'exact' grid 
results were obtained by repeated loops around the final V cycle of the FAS procedure, until 
the residuals were on the order of machine precision zero. The coarsest (full domain) scale 
had 5 points on a side, the next two finer scales 9 and 17, and the two nested patches both 
had 9 points on a side. To examine the order of the method, computations were performed on 
a full domain coarse grid corresponding to the that of the composite mesh (level h3), followed 
by one finer full domain grid with the spacing halved. The accuracy of the composite mesh 
method was then determined by comparison with the high order coarse and fine uniform 
grid calculations used in the determination of the order. 

The increased accuracy obtained using fourth order equations vs. second order is 
displayed for uniform domain computations in Figure 4, in which the absolute errors of the 
solution are presented away from the singularity. That the fourth order Laplacian leads 
to fourth order behavior was confirmed in the uniform grid computations described above. 
Except for the set of points one grid spacing away from the singularity, the correct order is 
obtained over the entire domain. 
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Then computations were performed on the four level composite domain with a single 
refinement patch centered at the origin. To test the effect of the boundary correction on 
conservation, the integral of the defect correction over the refinement patch was computed 
with and without the boundary terms. Without the boundary correction, the integral was 
0.8 in magnitude, while with the boundary terms, the integral was zero to double precision 
accuracy. The impact of the conservative boundary correction on the accuracy of the solution 
is apparent in Figure 5; serious errors are incurred over the whole domain in the absence of 
the boundary corrections. The accuracy of the method can be determined by comparison 
with the separate fine and coarse uniform domain results. Figure 6 shows that the accuracy 
within the patch is virtually identical to that for the uniform fine domain results. In Figure 7, 
the errors outside the refinement patch are displayed. The refinement mesh leads to increased 
accuracy outside the refinement on the coarser level in comparison with the uniform coarse 
level results. The numerical results are presented in Table 4. Finally, test computations 
were also successfully performed on a five level problem with two nested refinement patches. 
The resulting potential is plotted in Figure 8 to illustrate the accuracy of the method in 
relation to the numerical errors. These numerical results thus confirm that the conservative 
mesh refinement technique developed here leads to results of the correct high order within the 
refinement region, while increasing the accuracy on the coarse domain outside the refinement 
zone. 



6 Summary 

A general technique has been presented for carrying out high order mesh refinement multi- 
grid calculations. The FAS method for composite domains was first summarized. Then, 
Hamming's direct method for generating high order formulas was outlined. Both the high 
order Laplacian and interpolation coefficients were obtained from the universal matrices of 
the direct method. Since the sum of the defect correction over the interior of the patch 
is nonzero, high order conservative forms were derived by analysis of the one dimensional 
problem. The two and three dimensional forms can be obtained by averaging locally over 
three points on a line or nine points on a square, respectively. The new method was success- 
fully tested for the fourth order case on a Poisson problem, the source singularity in three 
dimensions. 

The high order mesh refinement methods should allow for accurate computations 
on three dimensional domains which require a range of length scales. We are develop- 
ing a quantum chemical Density Functional Theory (DFT) multigrid method for ab initio 



calculations. JT3[ So far, our fully numerical three dimensional calculations have been per- 
formed on uniform domains, treating both the electrons and nuclei with the high order grid 
approximations. As a test computation, we examined the CO molecule, and obtained good 
results in all electron computations. However, it is apparent from those results that the crude 
treatment of the nuclei and the core electrons limits the accuracy of the method. We plan to 
incorporate the high order composite mesh techniques into the quantum chemistry method 
to obtain more accurate results in the region of the nuclei, and to investigate the impact 
of those improvements on the eigenfunctions, eigenvalues, and total molecular energies. Be- 
yond the quantum chemical applications, these methods should prove helpful in other large 
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scale electrostatics calculations in biophysics and in multigrid fluid dynamics computations, 
especially for cases where increased resolution is required only over small local regions of 
space. 
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2nd 
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-490 


N=9 
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5040 


-9 


128 


-1008 


8064 


-14350 



Table 1: Coefficients for the Laplacian. One side plus the central point are shown. Each 
coefficient term should be divided by the prefactor. The Laplacian is symmetric about the 
central point. 
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Table 2: Coefficients for interpolation. One side of the symmetric weight vector is shown. 
Each coefficient term should be divided by the prefactor. 
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-9 


110 


-770 


6286 



Table 3: Coefficients for conservative forms. One side is shown. Each term on the other half 
of the gradient operator has the opposite sign. The set of coefficients is for a left boundary. 
All the signs are reversed for a right boundary. The locations of the terms for a 6th order 
example are shown in Fig. 3. 
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Exact 


Grid Exact 


Trunc. Err. 


Fine Trunc. Err. 


Coarse Trunc. Err. 


MG Err. 


2 


0.5 


0.520842737 


0.020842737 


0.020888618 




0.020875723 


4 


0.25 


0.257687505 


0.007687505 


0.007710007 


0.010443407 


0.007726450 


6 


0.16 


0.167993560 


0.001326893 


0.001183652 




0.001389060 
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0.125 


0.126091267 


0.001091267 


0.000232278 


0.003854094 


0.001232428 


12 


0.083 


0.0834252199 


0.0000918866 


0.0000295553 


0.0005908833 


0.0001265184 


16 


0.0625 


0.0625022779 


0.0000022779 


0.0000072359 


0.0001151041 


0.0000094832 


20 


0.05 


0.0499998476 


-0.0000001524 


0.0000023526 


0.0000347405 


0.0000016978 


24 


0.0416 


0.0416667925 


0.0000001258 


0.0000008825 


0.0000131217 


0.0000007387 



Table 4: Numerical results for the FAS-MG composite mesh computations at several dis- 
tances from the singular source. The edge of the patch is at r = 8. The last column is for a 
single V cycle MG computation with a total of 16 relaxation sweeps on the fine scale. Notice 
that for the points in the outer regions of the domain, the single cycle MG errors are not 
strictly less than the anomalously small truncation errors for the composite domain; however, 
the errors are considerably smaller than the truncation errors on the uniform coarse scale 
domain (they are of the same magnitude as the uniform fine grid results at those points). 



Figure 1: Schematic two dimensional cut through the three dimensional composite mesh. 
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Figure 2: Four level FAS cycle. 
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Figure 3: Locations and values for the coefficients used to generate the 6th order conservative 
form on a left boundary. The boundary is located on the coarse scale H at the point % and 
on the fine scale h at the point j = 2i. 
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Figure 4: Uniform domain calculations illustrating the improved accuracy of the fourth 
order approximation over the second order case. 
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Figure 5: Impact of nonconservation on the accuracy. The refinement patch edge is at r = 8. 
The diamonds are for the nonconservative calculation and the crosses for the conservative 
case. 
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Figure 6: Comparison of the accuracy for the composite grid with that on the uniform fine 
and coarse scales. The refinement patch edge is at r = 8. The diamonds are for the composite 
grid, the crosses for the uniform fine grid, and the squares for the uniform coarse grid. Note 
the potential value at the r = 8 point on the composite grid is set by the coarse grid, which 
leads to the larger error there. The interior patch points are at r = 2, 4, 6. 
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Figure 7: Comparison of accuracy outside the refinement patch with the uniform coarse 
grid results. The crosses are the uniform coarse grid results and the diamonds are for the 
composite grid computation. 




Figure 8: Plotted are the analytical 1/r potential and the numerical results from the con- 
servative mesh refinement multigrid computation. The two fine patches span the ranges -8. 
to 8. and -4. to 4. The lower curve gives the magnitude of the difference between the exact 
and numerical results, illustrating the larger errors near the source singularity. 
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